load ./libgeoconstraint.so.0.0.0
package require geoconstraint
proc calcforces { } {
    global atoms numatoms masses  namd  namd COM forceconstant groups 
    #load coorinates
    loadcoords coords
    #Determine center of mass from each group of atoms
    #groups of atom index orginized like that { { atom1 atom2 atom3} group2 } 
    foreach group $groups {
        set comsum {0.0 0.0 0.0}
        set totalmass 0.0
        foreach atom $group {
        set tmp [vecscale $masses($atom) $coords($atom)]
        set comsum [vecadd $comsum $tmp]
        set totalmass [expr $totalmass +$masses($atom)]
        }
        set com [vecscale [expr 1.0/$totalmass] $comsum]
        lassign $com x0 y0 z0
        # load Constraint value pre-calulated grids.
        set tmp [getConstraint geo $forceconstant $x0 $y0 $z0 v]
        set ener [lindex $tmp 0]
        set force [lrange $tmp 1 end]
        foreach atom $group {
            #separate the force to individual atoms
            set iforce [vecscale [expr $masses($atom)/$totalmass)] $force]
            addforce $atom $iforce 
        }
    }
}
